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Abstract 

Aspects of the QCD parton densities are briefly reviewed, drawing some paral- 
lels to the density matrix formulation of quantum mechanics, exemplified by Wigner 
functions. We elaborate on the solution of their evolution equations using loga- 
rithmic expansions and overview their kinetic interpretation. We illustrate how a 
Fokker-Planck equation can be derived using the master formulation of the same 
equations and its construction in the case of the transverse spin distributions. A 
simple connection of the leading order DGLAP equation to fractional diffusion using 
fractional calculus is also briefly outlined. ^ 



^To appear in Lecture Notes of SIM, Seminario Interdisciplinare di Matematica, S. Dragomir 
(ed.) 
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1 Introduction 



In Quantum Chromodynamics, or QCD, the accepted theory of the strong interactions, 
the description of high energy coHisions is formulated with the help of some important 
theorems, commonly refered to as factorization theorems According to these theorems 
it is possible to describe certain hadronic processes in three steps, the first involving 
the initial state, followed by a second intermediate state, commonly known as "the hard 
scattering", and the third which involves the final state. Both the initial and the final state 
are described with the help of some non-local operators supported on the light cone, whose 
matrix elements, termed "parton distributions" and "fragmentation functions" respectively, 
can't be computed from first principle, at least at this time, in a second quantized theory. 
It is possible, however, to parametrize their functional form by combining a vast amount 
of experimental data with numerical fits. 

The separation between the initial and the intermediate state is perfomed with the 
introduction of an intermediate scale, the factorization scale, denoted by Q, where Q is 
the amount of energy/virtuality involved in the actual separation. 

Parton distributions depend both on the choice of Q and on a second variable, x 
(Bjorken x) which is supported on the interval (0, 1). This variable can be interpreted as 
the fraction of the momentum of the particles in the initial state that fiows into the inter- 
mediate state. Parton distributions, at least some of them, have a simple interpretation in 
terms of densities of quarks and gluons (the partons) in a nucleon. They are probability 
distributions that change (evolve) with the scale Q and their evolution is described by a 
set of equations known as DGLAP equations ^, which are a special type of renormaliza- 
tion group equation (RGE). They are, by now, textbook material in high energy physics. 
Mathematicians may find a thorough discussion of their origin in [ij and in other books 
in Quantum Field Theory and particle phenomenology. 

On the contrary of other RGB's, such as those describing the evolution of other local 
operators in the Standard Model which are ordinary differential equations, the RGB's of 
the parton densities are integro-differential equations. They are defined in terms of some 
kernels {P{x)) which are not regular functions but distributions, finite after integration 
over X. 

The objective of this work is to present a short review of previous studies by us on this 
subject. Here we will discuss first a direct way to solve these equations and then move 
to describe their kinetic interpretation. We will show how one can construct a hierarchy 
of equations starting from the original integro-differential equations and illustrate the 
procedure at the lowest order and identify a Fokker-Planck approximation to the same 
equations. A simple connection between the leading order equations and the formalism 
of fractional calculus is also briefly pointed out. 

2 Parton distributions as Wigner functions 

Bvolution equations describing the high energy behaviour of scattering amplitudes carry 
significant information on the factorization/ renormalization scale dependence of such 
amplitudes, and allow to link the behaviour of processes at a given energy scale to collisions 
taking place at another (usually much higher) scale. 

In QCD confinement forbids the detection of the fundamental states of the theory, 
such as quarks and gluons. However, asymptotic freedom, the smallness of the coupling 
constant as the energy raises (or approaching smaller distances), allows to separate the 

^Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 
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perturbative dynamics at short distances from the non-perturbative one, due to confine- 
ment, through factorization theorems and the introduction of the parton distributions. We 
will very briefly introduce them below, and we will exploit the density matrix formulation 
of quantum mechanics as an analogy to illustrate the topic. 

The mathematical construct which is the closest to a parton distribution function 
(p.d.f.) q{x,Q) is a Wigner function, as first noticed in PJ. The analogy is, of course, 
limited. 

We recall that Wigner's description of quantum mechanics via quasi-probabilities of 
phase space (x, p) 

/(x, P) = ^ / dy^*(x - |y)e'P-''^(x + |y) (1) 

is fully equivalent to Schrodinger's formulation 0. Its dynamics is specified by Moyal's 
equation, that extends Liouville's formulation of classical mechanics for a classical hamil- 
tonian iJ(x, p), dtf + {/, H} = to quantum mechanics 

dl _ H*f-f*H 

dt ~ ilL ^ ' 



where the star-product (*) is 



A-i d ^dr,-d „d : 



^-\-^-v--^'- V. (3) 



The nonlocality of the construct appears immediately in its dynamics, with derivatives in 
its evolution equation Q that extend to all orders. In practice one can use translation of 
its arguments 

/(x, p) * g{x,p) = f (^x + z Y p, p - ^ Y 5^^^ 9{x, P) (4) 

for the explicit evaluation of this product. 

Differently from Wigner functions, in a p.d.f. the variable x takes now the role of 
the momentum "p". Also, parton distributions are correlation functions of a special type, 
being defined just on the light-cone. In this sense they are not generic nonlocal correlators. 
This limitation, due to the special nature of high energy collisions in asymptotically free 
theories, sets the boundary of validity of the parton model approach to QCD. 

Another observation that may help in the distinction between a p.d.f. and a Wigner 
function is the different nature of the respective equations. For the p.d.f. the nonlocality 
in the evolution is all expressed in a renormalization group equation rather than in the 
Hamiltonian formalism of However, also in this case some interesting analogies re- 
main. Once the evolution equation of a p.d.f. is written in the form of a master equation, 
a hierarchy of equations can be identified quite straighforwardly using the Kramers-Moyal 
expansion of the transition probabilities. This point will be discussed next. From this 
analysis, that we believe has limited phenomenological applications, emerge however some 
amusing features which may be generic for other types of nonlocal operators in field theory. 

2.1 The "Wiggies" of QCD 

In QCD one starts by introducing, via arguments based on unitarity, the hadronic tensor, 
which is the key construct describing the collision 
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^'^^ = / ^^"^^l [^''^o)' '^'(^)] 1^^*^^' 

with Pa and Pb being the momenta of the colliding hadrons and Sa and Sb their 
covariant spins. The J's are electromagnetic currents. 

The distribution functions that emerge -at leading order- from this factorized picture 
are correlation functions of non-local operators in configuration space. They are the 
quark-quark and the quark-antiquark correlators. 

Their expression simplifies in the axial gauge, in which the gauge link is removed by 
the gauge condition. For instance, the quark-quark correlator takes the form 

{<^^,a)^, (P, S,k) = j -0j-,e''-{P, S\4^\0)^|Ji^\z)\PS). (6) 

In (jni) we have included the quark flavour index a and an index A for the hadron, as 
usual. Fields are not time ordered since they can be described by the good light cone 
components oi and by At, a transverse component of the gauge field A^, as discussed 
in |4j. The non-perturbative information in a collision is carried by matrix elements of 
this type. 

Further considerations allow to show that the leading contributions to ^ come from 
the light-cone region. The leading expansion of the quark-quark correlator then is of the 
form 

I ^e''^PS\m^{Xn)\PS) = \ (Miix) + A75 hi{x) + 75 i^h,{x)) , (7) 

where we have used all the four vectors at our disposal (spin S, momentum P of the 
hadron) and introduced invariant amplitudes (parton distributions) /i, (71, /ii, now ex- 
pressed in terms of a scaling variable x (Bjorken x). is a light-cone four- vector (n^ = 0), 
approximately perpendicular to the hadron momentum. 

The definition of p.d.f.'s in ^ involves also an underlying physical scale Q {Q » 
^QCD, with Aqcd being the scale of confinement), not apparent from that equation and 
characterizing the energy scale at which these matrix elements, summarized by 0, are 
defined. Truly: /i = fi{x,Q^), hi = hi{x,Q^) and so on. 

The role of the various renormalization group equations is to describe the perturbative 
change in these functions as the scale Q is raised (lowered) . Each equation involves kernels 
(P(x)) of various types, of well known form, and asymptotic expansions of the solutions 
exist (see for instance [H] and the implementation given in jHI). The structure to all orders 
of the solution has been discussed by us recently using alternative methods 

2.2 Factorization and Evolution 

Factorization theorems play a crucial role in the application of perturbation theory to 
hadronic reactions. The proof of these theorems and the actual implementation of their 
implications has spanned a long time and has allowed to put the parton model under a 
stringent experimental test. Prior to embark on the discussion of our contributions to the 
study of evolution algorithms in Bjorken x-space, we provide here a brief background on 
the topic in order to make our treatment self-contained. We refer to [T] for a thorough 
overview of the topic. 

In sufficiently inclusive cross section, leading power factorization theorems allow to 
write down a hadronic cross section in terms of parton distributions and of some hard 
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scatterings, the latter being calculable at a given order in perturbation theory using the 
fundamental QCD lagrangean. Specifically, for a hadronic cross section, for instance a 
proton proton cross section app, the result of the calculation can be summarized by the 
formula 

^pp = ^ dxi I dx2fhi^f{xi,Q'^)fh2^f{x2,Q'^)a{xi,X2,s,i,Q'^) , (8) 

j: Jo Jo 

where the integral is a function of some variables Xi and X2 which describe the QCD 
dynamics at parton level in the Deep Inelastic Scattering (DIS) limit or, equivalently, 
at large energy and momentum transfers. These variables are termed "Bjorken variables" 
and are scale invariant. This formula is a statement about the computability of a hadronic 
collision in terms of some "building blocks" of easier definition. 

The variable Q^, in the equation above, can be identified with the factorization scale of 
the process, a can be computed at a given order in perturbation theory in an expansion in 
a;,,, the QCD coupling, while the fh.^f{x, Q"^) are the parton distributions. These describe 
the probability for a hadron h to prepare for the scattering a parton /, which undergoes 
the collision. An equivalent interpretation of the functions fh.^f{x,Q'^) is to characterize 
the density of partons of type / into a hadron of type h. A familiar notation, which 
simplifies the previous notations shown above, is to denote by qi{x, Q"^) the density of 
quarks in a hadron (a proton in this case) of flavour i and by g[x, Q^) the corresponding 
density of gluons. For instance, the annihilation channel of a quark with an antiquark in 
a generic process is accounted for by the contribution 



/ dxi 
Jo Jo 



dx2q{xi, Q'^)q{x2, Q'^)agg{xi, X2, Q^), (9) 



and so on for the other contributions, such as the quark-gluon sector (qg) or the gluon- 
gluon sector (gg) each of them characterized by a specific hard scattering cross section 
aqg, or aqg. In this separation of the cross section into contributions of hard scatterings a 
and parton distributions f{x, Q^) the scale at which the separation occurs is artificial, in 
the sense that the hadronic cross section a should not depend on or 

0. (10) 



dQ 



However, a perturbative computation performed by using the factorization formula, how- 
ever, shows that this is not the case, since the perturbative expansion of a 

a = a(°) + a,(QY^« + a,(gV^ (11) 

naturally brings in the dependence on the factorization scale Q. This dependence is weaker 
if we are able to push our computation of the hard scattering to a sufiiciently high order 
in ag. The order at which the perturbative expansion stops is also classified as a "leading 
order" (LO), "next-to-leading order" (NLO) or, even better, a "next-to-next-to-leading" 
(NNLO) contribution if more and more terms in the expansion are included. 

At the same time, as we have already mentioned, the parton distributions f{x, Q"^) 
also carry a similar dependence on the scale Q, which is summarized by some RGB's. The 
equations resum the logarithmic violations to the lowest order scale invariance, induced by 
the perturbative expansion. Also in this case we need to quantify this effect and reduce its 
impact on the prediction of the cross section. Solving the RGB's for the p.d.f.'s to higher 
order and, at the same time, computing the hard scatterings to higher orders reduces the 
spurious dependence on the factorization scale Q and improves the theoretical prediction 
of the real physical phenomenon. 
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2.3 Parton dynamics at NLO 



As we have mentioned, within the framework of the parton model, where the "partons" are 
the quarks and gluons inside a hadron, evolution equations of DGLAP-type - and the cor- 
responding initial conditions on the parton distributions - are among the most important 
blocks which characterize the description of the quark-gluon interaction. Other parts of 
this description require the computation of the hard scatterings (what we have called "the 
second stage" of the collision) with high accuracy. Here we illustrate an implementation 
to NLO of a method based on an ansatz which allows to rewrite the evolution equations 
as a set of recursion relations for some scale invariant functions, A„(x) and Bn{x), which 
appear in the expansion. The advantage, compared to others, of using these recursion 
relations is that just few iterates of these are necessary in order to obtain a stable so- 
lution. One of the advantages of the method is its analytical base, since the recursion 
relations can be solved explicitely in terms of the initial conditions. We have recently 
shown conclusively that these methods are equivalent to other methods that also give the 
exact solution using Mellin transforms. Our approach follows closely the Mellin method 
in the demonstration of the correctness of the recursion relations, but then is implemented 
directly in x-space. 



3 Running coupling 



In QCD the coupling constant changes with the change of the energy scale and there is a 
RGE associated to it as well, computed within a certain level of accuracy, also classified as 
LO, NLO, NNLO etc. The solution of the RGE for the running coupling constant in the 
NLO approximation can be expressed in terms of two coefficients (3q and describing 
the beta function of the theory and a renormalization group invariant scale Aqcd that 
sets the scale for confinement. The NNLO extension of the equation involves a third 
coefficient, ^2- The RGE of the coupling is 



da,(Q2 



d log ' 

and its expansion, for instance to NNLO, is given by [IH] [TT] 



Po 2 
Air ^ 



/?2 



IGvr 



2 s 



647r3 



at + 0(al 



(12) 



(13) 



For the moment we work in the NLO approximation (only (3o and /3i are included), and 
in this case an explicit solution of this equation is given by 



(3o \og{QVAl^) 



/3iloglog(QVA^ 



MS' 



Pi log(gvAk) 



+ 



log2(QVA2 



MS' 



(14) 



where 



and 



c 



Pi 



34 ,,2 
— A^, 



c 



^Ncnf-2CFnf, 



N'c 



4 
3' 



Tf = TrUj 



1 



(15) 



(16) 



2Nc 3' ' " ' 2 

where Nc is the number of colors, Uf is the number of active fiavors, which is fixed by the 
number of quarks with niq < Q. Aqcd is a scale that depends on the number of active 
quarks included in the evolution (number of fiavours) . We have taken for the quark masses 
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(charm, bottom and top respectively) rric = l.SGeV, = 4.5 GeV and rrit = 175 GeV, 
these are necessary in order to identify the thresholds at which the number of flavours n / 
is raised as we increase the final evolution scale. 

Aqcd should be more correctly denoted as A^^\ given its flavour dependence, and is 
given by 

^g4,5,6) ^ Q 248, 0.200, 0.131, 0.050 GeV. (17) 



The label "MS"' denotes conventionally the modified minimal subtraction scheme and 
denotes a specific regularization scheme of the perturbative expansion. 

We also define the distribution of a given helicity (±), /^(x, Q^), which is the prob- 
ability of finding a parton of type / at a scale Q, where / = qi,qi,g, in a longitudinally 
polarized proton with the spin aligned (+) or anti-aligned (-) respect to the proton spin 
and carrying a fraction x of the proton's momentum. 

We introduce the longitudinally polarized parton distribution of the proton 

A/(x,g2) = /+(x,g2)-r(x,g2). (is) 

We also introduce another type of parton density, termed transverse spin distribution, 
which is defined as the probability of finding a parton of type / in a transversely polarized 
proton with its spin parallel (|) minus the probability of finding it antiparallel (|) to the 
proton spin 

ATf{x,Q')^f\x,Q')-f^{x,Q'). (19) 
Similarly, the unpolarized (spin averaged) parton distribution of the proton is given 

by 

fix, Q') ^ f^ix, Q') + fix, Q') = fix, Q') + fix, Q'). (20) 

We also recall, if not obvious, that taking linear combinations of Equations (f20|) and (fT8| . 
one recovers the parton distributions of a given helicity 

;.(,.Q.),/(^.Q')±A/(..Q-)^ (21) 

In regard to the kernels, the notations P, AP, A^P, P"*"^, will be used to denote the ker- 
nels in the unpolarized, longitudinally polarized, transversely polarized, and the positive 
(negative) helicity cases respectively. 

The DGLAP equation is an integro-differential equation whose general mathematical 
structure is 

fix, Q') = Pix, asiQ')) ® fix, Q'), (22) 



dlogQ2 

where the convolution product is defined by 



® b] (^) = / -f^[-) b{y) = I -faiy)b ^- ) . (23) 



mg 



y \yj Jx y 

Let us now turn to the evolution equations, starting from the unpolarized case. Defin- 

^ =g^±q^, q^^^ = J2 ^^^^ ' = ^^^^ - -^^^^ ' (24) 



i=l 

the evolution equations are 

d 



dlogQ 



\-\x,Q^) = PNS-{x,asiQ^))®q^\x,Q^), (25) 
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^ Xrix,Q') = PnsA^,c^s{Q')) ® Xi{x,Q'), (26) 



d log 

for the non-singlet sector and 

d /gW(x,g^)\_/P,,(x,a,(g^)) n,(x,a,(Q2)) \ /g{+)(^,g2) X 
dlogQ2^ 9ix,Q') J [Pg,{x,asm) P^^^asm)) [ 9ix,Q') J 

for the singlet sector. Equations analogous to ()24ll27p . with just a change of notation, 
are valid in the longitudinally polarized case and, due to the linearity of the evolution 
equations, also for the distributions in the helicity basis. In the transverse case instead, 
there is no coupling between gluons and quarks, so the singlet sector (|27|l is missing. In 
this case we will solve just the nonsinglet equations 

;ATqt\x,Q') = ATPNsAx,asm) ® ATqt\x,Q^), (28) 



diogg2 

^ ;ATqt\x,Q') = ATPNS+{x,as{Q')) ® ATqt\x,Q'). (29) 



dlogQ2 

We also recall that the perturbative expansion, up to next-to-leading order, of the 
kernels is 

P{x,a.) = (g) P^'\x) + {l^yP''\x) + .... (30) 

Kernels of fixed helicity can also be introduced with P++{z) = {P{z) + AP{z))/2 and 
P+_(^) = {P{z) — AP{z))/2 denoting splitting functions of fixed helicity, which will be 
used below. 

The equations, in the helicity basis, are 

at 2n y y 

+PIM-) ® g+iy) + Pl'-{-) ® 9-iy)), 
y y 

= ^iP^-i-) ® + n+(-) ® <i-{y) 

at 2Ti y y 

+prC-)®9+{y) + Pl%C-)®9-{y)). 
y y 



+pf+(-) ® 9M + Pi'-i-) ® 9-iy)), 
y y 

+Pf _(-) ® gM + PIM-) ® 9-{y)). (31) 



The non-singlet (valence) analogue of this equation is also easy to write down 

= ^(^++(-) ® Q^yiy) + P^-i-) ® Q-yiy)), 

at 271 y y 

= ^(n-(-) ® <i+,v{y) + ® <i-,v{y))- (32) 

at 2tt y y 
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where the q±y = q± — q± are the valence components of fixed helicity. The kernels in this 



basis are given by 






p(0) 

^NS±,++ 


_ p(0) 


+ ~ ^qq 


p(0) 

-^qq,+- 


_ p(0) 


+ = 


p(0) 

q9,++ 


= Ufx'^ 




Pq9,+- 


— P 

qg- 


-+ = nf{x 


Pgq,++ 


= p 
^ gq- 


1 

. = C,^ 

X 


p(0) 


_ p(0) 

ggA 


/ 


p(0) 


= nJ 


3x + - - 

X 



I? 



^1 — x)+ X 

o ™2 



+ --l-x-x^] +(3o6{l 



(33) 

Taking linear combinations of these equations (adding and subtracting) , one recovers 
the usual evolutions for unpolarized q{x) and longitudinally polarized Aq{x) distributions. 



4 The mathematical structure of the kernel 

Here we try to illustrate some simple manipulations on the kernels which are very useful 
in order to simplify the equations. Their typical form is the following 

Pix) = P,{x) + + P,5{1 - x), (34) 

with a regular part P\{x^, a "plus distribution" part Piix^ji^- — a;)+ and a part P3 which 
multiplies a Dirac delta function. For a generic function a(a;) defined in the [0, 1) interval 
and singular in a; = 1, the plus distribution [a(x)]+ is defined by 

f{x)[a{x)]+dx= [ {f{x)-f{l))a{x)dx, (35) 

where f{x) is a regular test function. Alternatively, an operative definition (in the sense 
of distributions) is the following 

[a(x)]_|_ = a{x) — 6{1 — x) a{y)dy. (36) 







From it follows immediately that each plus distribution integrate to zero in the [0, 1] 
interval 

• 1 

[a{x)]+dx = 0. (37) 



We want to make the convolution of the generic kernel ([T?!!) with a function f{x). 
We introduce the notation /(x) = xf{x), the factor x being introduced to stabilize the 
evolution for small x values. The treatment of the regular and the delta-function parts is 
trivial 

Pi(x) ® fix) = xP,{x) ® fix) =x f ^Pi{y)f = £ dyP^iy)f (^^^ (38) 



X 



Jx y \yj 



P,S{1 -x)® fix) = / -^PsSil - y)/ - . (39) 
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Let us now treat the more involved case of the plus distribution part 



(1 — x)+ 1 — X 



^-yJx y \yj 



f - 

z y ^-y \y 
y ^-y \yj Jo i -y 

2/1-?/ \yj Jx 1-2/ 

-A(l)/(x)£^ (40) 

^ d.p.(2/)/(x/2/)-2/p.(l)/(x)°^ ^^^^ _ (^^^ 



which yields 



P2(y)/>/2/)-P2(l)/» , , 

®f{x)= dy + /(a;)log(l - x). (42) 



(1 l-y 

From Feynman diagrams calculations one can get just the regular part Pi{x) of each 
kernel. The remaining distributional parts (plus distribution and delta distribution) 
emerge from a procedure of regularization, that introduce the plus distribution part to reg- 
ularize the eventual singularity in a; = 1 and the delta distribution to fulfil some physical 
constraints, the sum rules. 

The first one is the baryon number sum rule (BNSR), asserting that the baryon number 
(number of quarks less number of antiquarks) of the hadron must remain equal to its initial 
value (3 in the case of the proton) throughout the evolution, i.e. for each value of 



Jo 



q^->{x,Q^)dx = 3. (43) 

Deriving with respect to log(Q^) and having in mind that q^~^ evolves with P^^, we 
get 

[\x[P);s{Q')®q^-\Q')]ix) = 0. (44) 
Making use of the property of the Mellin moment of a convolution (TTHIl this implies 

(^^' P^5(x, Q')dx^ (^^' g(-)(x, Q')dx^ = 0, (45) 

from which, using (|i3|) . we find the BNSR condition on the kernel 



P^5(a;)da; = 0. (46) 







The other constraint is the momentum sum rule (MSR), asserting that the total mo- 
mentum of the hadron is constant throughout the evolution. Having in mind that x is 
the fraction of momentum carried out by each parton, this concept is translated into the 
relation ^ 

[ (xg(+)(x, Q^) + xg{x, Q^)) dx = 1 (47) 
Jo 
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that must hold for each value of . Deriving with respect to log((5^) and using the 
singlet DGLAP equation one obtains 



In a similar fashion, using (fTBl we also obtain 



(48) 



+ 



{P,g{X,Q^)+Pgg{X,Q^))dX 



{Pgg{X,Q^)+Pgg{X,Q^))dX 



xq^~^\x, Q^)dx 



xg{x, Q^)dx 



(49) 



from which we find the MSR conditions on the singlet kernels 



X 



X 



{P,g{X,Q^) + Pg,{X,Q^)) dX = 0, 
{Pgg{X,Q^)+PggiX,Q^)) dX = . 



(50) 
(51) 



4.1 The kernels and their regularization 

We illustrate now an example of the regularization procedure of the DGLAP kernels 
through the sum rules. The LO kernels computed by diagrammatic techniques for x < 1 
are 



PnUx) = Cp 



l+x^ 



pS\x) = 2N, 



1 — xj "[1 — X 

Pi^{x) = 2Tf[x' + {l-xr] 

pr(x)=c/'+^'-"'' 
1 



1 — X 



X 



c 



1 — X X 



H 2 + x(l - x) 



(52) 
(53) 
(54) 

(55) 



We want to analytically continue these kernels to x = 1 curing the ultraviolet singularities 
in Pi? (a:) 

and Pgg\x). We start introducing the plus distribution prescription in Pqq\x). 
We make the replacement 

11 

1^ ~" (l-x)+ ^^^^ 

to avoid the singularity and we add a term M(l — x) (where k has to be determined) to 
fulfill the BNSR So we have 



X 



;i-x) 



- 1 - X + k5{l - x) 



+ 



(57) 



Imposing by the BNSR that Pq^{x) integrates to zero in [0, 1] and remembering that the 
plus distribution integrates to zero we get 



p(°)(x)dx = C^ 



-l-- + k 
2 



0, 



(58) 
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hence k = 3/2, and the regularized form of the kernel is 



- — 1-X + -6(1 - x) 

(l-x)+ 2 ^ ^ 



(59) 



Noticing that 



X 



Jo (1-3;)+^'^ JO 



-1 + 



dx = —1 



(60) 



it can be easily proved that the MSR (f5n|l is satisfied. Let us now regularize Pgg{x). We 
make the replacement 



c 



1 



1 



, 1 — X 



-2 + x(l-x) 



+ k6{l-x). 



(61) 



Imposing the other MSR (f5T| we get 

»i 



2N, 



c 



X 



[1 - x) + 



1 - 2x + x^(l - x) 



+ kx6{l — x) 



+2Tf [x^ + x{l - xY]} dx = 0, 



from which we find 



c - i^Tf 



f3o 



so the regularized form of the kernel is 

P^'\x) = 2Nc 



11 . ^ 

— + --2 + X 1-x 

1 — x)+ X 



+ ^S{l-x). 



(62) 
(63) 

(64) 



5 First view: an ansatz from x-space and some exam- 
ples 

In order to solve the evolution equations directly in x-space, we assume solutions of the 
form 



as(Q, 



2\' 

o) 



(65) 



n=0 „=o 

for each parton distribution /, where Qo defines the initial evolution scale. The justifica- 
tion of this ansatz can be found, at least in the case of the photon structure function, in 
the original work of Rossi ^2)1, and its connection to the ordinary solutions of the DGLAP 
equations is most easily worked out by taking moments of the scale invariant coefficient 
functions A„ and and comparing them to the corresponding moments of the parton 
distributions, as we are going to illustrate in section 5. The link between Rossi's expansion 
and the solution of the evolution equations (which are ordinary differential equations) in 
the space of the moments up to NLO will be discussed in that section, from which it will 
be clear that Rossi's ansatz involves a resummation of the ordinary Mellin moments of 
the parton distributions. 

Setting Q = Qo in (j65|) we get 



f{x,Ql) = Ao{x) + as{Ql)Bo{x). 



(66) 



12 



Inserting (fHK| in the evolution equations, we obtain the following recursion relations for 
the coefficients An and Bn 

An+i{x) = -^P^'\x) ® An{x), (67) 
Po 



= - -/^An+iix) - ^P^'Hx) ® Bn{x) - -^P^'\x) ® An{x) (68) 

47rPo Po vrpo 

obtained by equating left-hand sides and right-hand-side of the equation of the same 
logarithmic power in log"as((5^) ^ind log"^ ^^(Q^). Any boundary condition satisfying 
(|66|) can be chosen at the lowest scale Qo and in our case we choose 

Bo{x) = 0, f{x,Ql)=Ao{x). (69) 

In the numerical implementations of this procedure one has to be careful in the manip- 
ulations of the numerical integrals, being all the kernels defined as distributions. Since 
the distributions are integrated, there are various ways to render the integrals finite, as 
discussed in the previous literature on the method |T3| in the case of the photon struc- 
ture function. In these previous studies the edge-point contributions - i.e. the terms 
which multiply 6{1 — x) in the kernels - are approximated using a sequence of functions 
converging to the 6 function in a distributional sense. 

This technique is not very efficient. We think that the best way to proceed is to 
actually perform the integrals explicitly in the recursion relations and let the subtracting 
terms appear under the same integral together with the bulk contributions (x < 1). This 
procedure is best exemplified by the integral relation 

dy fi I \ r dyyf{y) - xf{x) 

-f{x/y) = / log(l - x)f{x) (70) 



y(i-y)+ Jx y y-x 

in which, on the right hand side, regularity of both the first and the second term is explicit. 
For instance, the evolution equations become (prior to separation between singlet and 
non-singlet sectors) in the unpolarized case 

dqi{x) f dy yqi{y) - xqi{x) , , ^ , ^ ^y 3 

- 2Cf / h 2Cirlog(l - x)gi(x) - / — [l + z) qi[y) + -C f 



d\og{Q^) J y y - X y 2 

"dy 



+nJ ^{z' + {l-zr)giy) 

J X y 



dg{x) r dy l + jl-zf dy yf{y)-xf{x) 

Cp qi{y) + 2Nc g{y) 



diog(Q2) y z Jx y y-x 

log(l - x)g{x) + 2iV, ^ ^ (^^ - 2 + z{l - z)^ g{y) + ^g{x) 



with z = x/y. Here q are fixed fiavour distributions. 



6 The evolution of the transverse spin distributions 

We show here an example, worked out in some detail, that illustrates the kind of manip- 
ulations that are involved in the computation of the recursion relations. 
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In fact, leading order (LO) and NLO recursion relations for the coefficients of the 
expansion can be worked out quite easily. We illustrate here in detail the implementation 
of a non-singlet evolutions, such as those involving transverse spin distributions. For the 
first recursion relation (f67|) in this case we have 



f3o 



^ dyyA^{y)-xAi 



X y 



y-x 



+ A^{x) log(l-x) 



+ 



(72) 



As we move to NLO, it is convenient to summarize the structure of the transverse kernel 
ATPi'^^\x) [H] as 



/\TPt/^\x) = Ktix)5il -x)+ K^{x)S2{x) + Kf{x) log(x) 



+K^{x) log^(x) + K^{x) \og{x) log(l -x) + K^{x) 



Hence, for the (+) case we have 



+ Kf{x). (73) 



(1-x) 



J X y 

+ log^z)K+(z) + log(z) log(l - z)K+(z)] Ai(y) + 

vj X y y X ) Jx y 

where z = x/y. For the (— ) case we get a similar expression. 
For the B^j^^{x) we get (for the (+) case) 



B+^,{x) = -Bt{x) + ^,\2CF 



^ dy yA+{y)-xA+{x) 

y y-x 



+ A+(x)log(l-x) 



+ 



-2Cp jK{y)) + C/-Al{x)^ - ^KtAl{x) + £ ^ [Kt{z)S,{z) + 



+ Kt{z) \og{z) + \og\z)Kt{z) + log(z) log(l - z)Kt{z 



1 



K{y) + 



^ dyyA+{:y)-xA+{x) 



4 

Jo 



X y y-x 

UyyBt{y)-xBt{x) 



+ A+(x)log(l -x) 



47r/?o 

+ Kt ^-^A^iy) 
J X y 



X y 



y-x 



+ B^{x) log(l-x) 



+ 



.Poj \Jx y 

The explicit expressions of the Ki can be found in [TK] . 
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6.1 Solutions in moments space and comparisons 

It is interesting to compare the recursive solution with the solution in moment space, that 
is easy to derive. In moment space the equations become algebraic and can be readily 
solved. For this we need some definitions. 

The n-th Mellin moment of a function of the Bjorken variable f{x) is defined by 

fn= f x''-^f{x)dx. (75) 

An important property of Mellin moments is that the Mellin moment of the convolution 
of two functions is equal to the product of the individual Mellin moments 

[f®9]n = fngn. (76) 

Let us prove it. 

[/ ® g]n = I' dx X-'' [f ® g] {x) = dx X--' £ ^f{y)g (^^^ . (77) 
Exchanging the x and y integrations 

[/ ® 9l = f f{y) r -x'^-'g (-) , (78) 
Jo Jo y \y / 

and introducing the new variable z = x/y 



1 p1 



[f®9l= / dyf{y)y''~' dz z--'g{z) = f^g^,. (79) 
Jo Jo 

This leads to an alternative formulation of DGLAP equation, that is also the most 
widely used to solve numerically the evolution equations. By taking the first Mellin mo- 
ment of both sides of the integro-differential equation (f22|l we are left with the differential 
equation 

fi{Q') = PimMQ') (80) 



dlogQ2 

that can be easily solved to give 



/i(Q')= / fix,Q')dx. (81) 





To get the desired solution /(x, Q^) there is a last step, the inverse Mellin transform of 
the first moment of the parton distributions, involving a numerical integration on the 
complex plane. This is the most difficult (and time-consuming) task that the algorithms 
of solution of DGLAP equation based on Mellin transformation - by far the most widely 
used - must accomplish. It is particularly instructing to illustrate here briefiy the relation 
between the Mellin moments of the parton distributions, which evolve with standard 
ordinary differential equations, and those of the arbitrary coefficient A„(x) and Bn{x) 
which characterize Rossi's expansion up to next-to leading order (NLO). This relation, 
as we are going to show, involves a resummation of the ordinary moments of the parton 
distributions. 
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Specifically, here we will be dealing with the relation between the Mellin moments of 
the coefficients appearing in the expansion 



1 

N-l 



A{N) = / dxx^-^A{x) 
Jo 

Jo 



(82) 



and those of the distributions 



ATq^^\N,Q^)= [' dxx''-'ATq^^\x,Q')). (83) 
Jo 

For this purpose we recall that the general (non-singlet) solution to NLO for the latter 
moments is given by 



ATq±{N,Q') = K{Ql,Q\N)\^^-^j ATq±{N,Ql) 
with the input distributions Axq^iQl) the input scale Qq. We also have set 
K{Ql Q\ N) = l+ "-(^o)^-^c..(g^) _ ^ArP^^^^N) 



2Po 



(84) 



In the expressions above we have introduced the corresponding moments for the LO and 
NLO kernels (AtP,?^'"^, ArP^l^f). 

The relation between the moments of the coefficients of the non-singlet x— space ex- 
pansion and those of the parton distributions at any Q, as expressed by eq. (|8Hl can be 
easily written down 

(—2 A P (N)\^ 
V ) • 

As a check of this expression, notice that the initial condition is easily obtained from 
setting Q ^ Qo,n ^ 0, thereby obtaining 

A^'{N) + a,B^'{N) = Arq^iN, Ql), (86) 

which can be solved with A^^N) = ATq±{N,Ql) and B^'^iN) = 0. 

It is then evident that the expansion (|65|) involves a resummation of the logarithmic 
contributions, as shown in eq. (|^ . 

In the singlet sector we can work out a similar relation both to LO 



,„,^H,|^)V,,(^, ,37, 



with 



ei = ^-^(P(°HiV)-A,l) 

A2 — M 



Ai,2 = I (^p!^^\N) + pW(iV) ± y [pS\n) - Pi?(iV))' + 4p(?(iV)pW(iV) j (88) 
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and to NLO 



where 



with 



A^{N) + a,B^{N) = xi { -^J + X2 (^-^ ] , (89) 



« /-2 e2Rei 
Xi = ei + — — eiRei + 



27r V /?o Ai - A2 - /3o/2 

a /-2 eiRes 
X2 = 62 + — — e2Re2 + 



2vr V /?o A2 - Ai - /3o/2 



(90) 



R = p(^\N) - Ap{o)(jv). (91) 
2po 

We remark, if not obvious, that An{N) and S„(A^), P(°)(A^), P^^'^N), in this case, are all 
2-by-2 singlet matrices. 



7 Moving to higher orders 

The structure of the solution to higher orders can be worked out in generality. More 
details can be found in 0, here we just outline the procedure. 

The perturbative expansion of the kernels now includes the NNLO contributions and 
is given by 

P{x,a.) = (^) P^'\x) + {^yP^'\x) + (l^yP^'Hx) + .... (92) 

whose specific form can be found in the original literature [T6| fT7]. 

We solve Eq. (f^ directly in x-space, assuming a solution of the form 



n=0 



for each parton distribution /, where Qq defines the initial evolution scale. 

As in the previous examples, also in this case we derive the following recursion relations 



for the coefficients An, Bn and C„ 



Bn+i(x) = -Bjx) - 



An+i{x) = -^P^'\x) ® An{x), (94) 
Po 

VA„+i(x) - ^P^'\x) ® Bn{x) - -^P^'\x) ® A„(x), (95) 



4vr/3o Po tx(3q 



Cn-,.{X) = -2C^n(^)-if^^n(^)-if^5n-M(^)-Yj^^n-.l(x) 

_ 2 p(o)^^) ® Cn{x) - -^P^'\x) ® 5„(x) 

-7ATP^'\x)^An{x). (96) 
27r^Po 
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It is an easy exercise to derive the recursion relations for the coefficients of the expansion. 
We illustrate the derivation for the interested reader. 
We introduce the notation 



(9^ 



and, making use of the beta function definition (fT2|l . we compute its derivative 



Inserting our ansatz ([HHll for the solution into the DGLAP equation we get for the 
LHS 



n=l 



n! Us n\ as 



OO f 

n=0 



+ 2a,/5(a,)^L" |. . (99) 
n\ TV. ' 



Sending n — > n — 1 in the first sum, using the three-loop expansion of the beta function 
(fTr?|) and neglecting all the terms of order or more, the previous formula becomes 



n=0 



^ I n\ \ Arr " 



2 3 



n\ \ An 167r2 * 647r^ 



n! 



47r 167r2 " n\ \ An 



Using the expansion of the kernels (jHSI), we get for the RHS 

n=0 ■ 



^ [P(^) ® A„] (x) + $ [P(°) ® 



47r2 



[P(i) ® P„,] (x) + ^ [P(°) ® a] (a;) I . (101) 



Equating ()inn|) and (|inHl term by term and grouping the terms proportional respectively 
to Og, a1 and we get the three desired recursion relations (f9i|) . (|95|l and (|96|) . 
Setting Q = Qq m (jnTHl we get 

/(x, Ql) = Ao{x) + a,(Q2)5o(x) + {as{Q^)f Co{x). (102) 

We solve these relations using a boundary condition satisfying p()2l) . In our case we 
choose 

B^{x) = Co(x) = 0, /(x, Ql) = A,{x). (103) 
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The procedure studied in the previous section can be generalized and applied to obtain 
solutions that retain higher order logarithmic contributions in the NLO/NNLO singlet 
cases. The same procedure is also the one that has been implemented in all the existing 
codes for the singlet: one has to truncate the equation and then try to reach the exact 
solution by a sufficiently high number of iterates. These arguments can be extended to 
even higher orders. We have recently shown [Z| that the exact solution to all orders of 
the RGE's of the p.d.f.'s can be written down in close form by the introduction of a 
generalized logarithmic expansion that includes contributions of all powers of the strong 
coupling constant and are summarized by the expression 



k=0 

is determined perturbatively by the recursione relations extracted from the RGE. The 
analysis is quite involved and can be found in the original literature. 

8 Other solutions in x-space: the Laguerre expansion 

Another possible way to solve the evolution equations is to use an operatorial formalism 
and observe that the distributions can be expanded in a suitable orthonormal basis (La- 
guerre polynomials). This method has had limited applications due to the difficulty to 
handle numerically the small-x behaviour of the algorithm that constructs the solution. 
The implementation of the method is discussed in [0]. One starts by writing the expres- 
sion of the solution in operatorial form in terms of two singlet evolution operators E±(t, x) 
and initial conditions q±{x,t = 0) = q±{x) as 




(104) 



where the function 



oo 




(105) 



(106) 



whose solutions are given by 



ql \t,x) = E^^)®q\ > 
Xiit,x) = E(+) (g) Xi{x). 



(107) 



The singlet evolution for the matrix operator E{t, x) 




(108) 



dE 



P(g)E 



(109) 



is solved similarly as 




(110) 



(111) 
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This method, proposed in [.8j, requires an expansion of the splitting functions and of 
the parton distributions in the basis of the Laguerre polynomials 

fe=0 ^ ^ 

which satisfies the property of closure under a convolution 

Lniy) ® Lmiy) = Ln+m{y) " ^„+m+l(Z/)- (113) 

In order to improve the small-x behaviour of the algorithm, from now on, the evolution 
is applied to the modified kernel xP(x), which, for simplicity, is still denoted as in all the 
equations above, i.e. by P{x). At a second step, the < x < 1 interval is mapped into 
an infinite interval < y < oo by a change of variable x = e^^ and all the integrations 
are performed in this last interval. We start from the non-singlet case by defining the 
Laguerre expansion of the kernels and the corresponding (Laguerre) moments to lowest 
order 



n=0 

I'OO 

p(o) = / dye-yiMP^'^iy) (114) 

^0 



and to NLO 



R{y) = Y,RnLn{y)- (115) 

n=0 

One defines also the difference of moments 

n = R^-R^-l R-1 = 0. (116) 
A similar expansion is set up for the evolution operators E(t,y) 



n=0 
oo 

E{t,y) = Y,En{t)Lniy), (117) 



n=0 

where all the information on the t evolution is contained in the moments En{t). The 
solution to NLO is expressed as 

£»W = s,?'W-|-^^«^£<"W, (118) 
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where 

fe=0 

n 

E^:\t) = Y,T^.,Ef\t), (120) 

(121) 

and the coefficients An'' are determined recursively from the moments of the lowest order 
kernel P(°) 

4^=^^^ = E^^nW (^ = 0,l,2,...,n-l). (122) 

i=k 

In the singlet case one proceeds in a similar way. The solution is expressed in terms 
of a 2-by-2 matrix operator 

00 

£;(°)(t,y) = E^^H^)i^n(y). (123) 

n=0 

The solution (at leading order) is written down in terms of 2 projection matrices and one 
eigenvalue (A) of the P^^^ (matrix) kernel 

ei = ^P(°\ = ^ (-P(°) + Al) , (124) 

(125) 



where 



A = -(^Cf + ^n/Ts), (126) 



in the form 



= E^K'^ + ^i'^^'O- (127) 

fc=0 

The recursion relations which allow to build An and Bn are solved in two steps as 
follows. One solves first for two sets of matrices an^ and hn"* by the relations 

ai°) = 



n 



a 



n-l 



i=k 



= 



n-l 



i=l 
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which are used to construct the matrices An'^ and Bn^ 



5f = ei + (eiai") - (-ire^fei")) • (129) 



(130) 



These matrices are then input in the recursion relations 



n-1 



^(fc+1) _ Xq^j\(^) _|_ y^^p^^_.A'^ 

i=k 
n-1 

5^ = -Ae.5f + (131) 



i=k 



witn n > and A; = 0, 1, ...,n — 1, which generates the coefficients of the matrix-valued 
operator E^'^^ (i.e. the leading order solution). The NLO part of the evolution is obtained 
from 



E«(t,y) = Y.^n\t)Uy), (132) 



n=0 



with 



where 



E^:\t) = E^\t) - 2E^^\(t) + E^^l.it) (133) 



£:(i)(t) = I rfre-*"/2 J]E(°)(t - r)i?,-4°V)'^(^ - z - j - A;). (134) 

The expressions of E^'^'^ and i?^^^ are inserted into eq. ()118p thereby providing a complete 
NLO solution of the singlet sector. The implementation of this method to NLO is quite 
involved and provides an accuracy which is comparable to other methods of solutions. 



9 Second view: the master form of the evolution equa- 
tions and a hierarchy 

A second part of our brief review is dedicated to the discussion of the kinetic structure 
that underlines an evolution equation of DGLAP type. We recall that master equations 
are introduced in statistical mechanics in order to describe the evolution in time of some 
diffusive processes. One can take as an example an ensemble of random walkers moving 
on a given lattice. The density of walkers present in a certain point on the lattice can be 
described by a master equation and expressed in terms of transition probabilities. The 
master form of the DGLAP equation can be derived quite simply by performing some 
manipulations on the kernels, having separated the bulk contributions {x < 1) from the 
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Figure 1: The constrained random walk of the parton densities 



edge points (x = 1). Below we show how this is obtained and how one can contruct, by a 
formal Kramers-Moyal expansion of the transition probabilities, a hierarchy of PDE's of 
arbitrarily high orders. 

Let's start considering a generic 1-D master equation for transition probabilities w{x\x') 
which we interpret as the probability of making a transition to a point x given a starting 
point x' for a given physical system. The picture we have in mind is that of a gas of 
particles making collisions in 1-D and entering the interval (x, x + dx) with a probability 
w{x\x') per single transition, or leaving it with a transition probability w{x'\x). In general 
one writes down a master equation 

d f 

—f{x,T)= / dx' {w{x\x')f{x' ,t) — w{x'\x)f{x,T)) dx' . (135) 



Or 

describing the time r evolution of the density of the gas undergoing collisions or the motion 
of a many replicas of walkers of density /(x, r) jumping with a pre-assigned probability, 
according to taste. 

The result of Collins a Qiu, who were after a derivation of the DGLAP equation 
that could include automatically also the "edge point" contributions (or x=l terms of the 
DGLAP kernels) is in pointing out the existence of a probabilistic picture of the DGLAP 
dynamics. These edge point terms had been always introduced in the past only by hand 
and serve to enforce the baryon number sum rule and the momentum sum rule as Q, the 
momentum scale, varies. 

The kinetic interpretation was used in [TH] to provide an alternative proof of Soffer's 
inequality. We recall that this inequality 

|/ii(x)| < g+(x) (136) 

famous by now, sets a bound on the transverse spin distribution hi{x) in terms of the 
components of the positive helicity component of the quarks, for a given flavour. The 
inequality has to be respected by the evolution. We recall that hi, also denoted by the 
symbol 

Arqix, Q2) = qHx, Q^) - q\x, Q^), (137) 

has the property of being purely non-singlet and of appearing at leading twist. It is 
identifiable in transversely polarized hadron-hadron collisions and not in Deep Inelastic 
Scattering (DIS), where can appear only through an insertion of the electron mass in the 
unitarity graph of DIS. 

The connection between the Collins-Qiu form of the DGLAP equation and the master 
equation is established as follows. The DGLAP equation, in its original formulation is 
generically written as 

(138) 

rflog(g^) y 

where we are assuming a scalar form of the equation, such as in the non-singlet sector. 
The generalization to the singlet sector of the arguments given below is, of course, quite 
straightforward. To arrive at a probabilistic picture of the equation we start reinterpreting 
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r = log((5^) as a time variable, while the parton density q{x, r) lives in a one dimensional 
(Bjorken) x space. 

We recall that the kernels are defined as "plus" distributions. Conservation of baryon 
number, for instance, is enforced by the addition of edge-point contributions proportional 
to (5(1 - x). 

We start with the following form of the kernel 

P{z)=P{z)-6{l-z) [ P{z)dz, (139) 

Jo 

where we have separated the edge point contributions from the rest of the kernel, here 
called P{z). This manipulation is understood in all the equations that follow. The 
equation is rewritten in the following form 



d_ 

dr 

Now, if we define 



qix, r) = C dyP fdj^p (I) ^ (140) 

Jx \yj y Jo y ^^J x 

wix\y) = -±P{x/y)^ ^ (141) 

ivr y 



()14n|) becomes a master equation for the probability function q{x, r) 
d f 

—q{x,T)= / dx' {w{x\x')q{x' ,t) — w{x'\x)q{x,T)) dx' . 



(142) 



There are some interesting features of this special master equation. Differently from other 
master equations, where transitions are allowed from a given x both toward y > x and 
y < X, in this case, transitions toward x take place only from values y > x and leave the 
momentum cell {x,x + dx) only toward smaller y values (see Fig.JH). 

Clearly, this sets a direction of the kinetic evolution of the densities from large x values 
toward smaller-x values as r, the fictitious "time" variable, increases. 

Probably this is the simplest illustration of the fact that parton densities, at large 
final evolution scales, are dominated by their small-x behaviour. As the "randomly moving 
partons" reach the x ~ region of momentum space, they find no space where to go, while 
other partons tend to pile up toward the same region from above. This is the picture of a 
random walk biased to move downward (toward the small-x region) and is illustrated in 
Fig. dH). 



10 Probabilistic Kernels 

We briefly discuss some salient features of the structure of the kernels in this approach 
and comment on the type of regularization involved in order to deflne them appropriately. 
We recall that unpolarized and polarized kernels, in leading order, are given by 

PSI = Pi^=CF{j^^-l-X+Pil-x)^ 

i^J) = 2Tj {x' + {1 - x^)) 
p(0) ^ ^,l±(i^ 



gg 



X 



= 2iV.(^^^ + i-2 + a:(l-a:))+|5(l-x) (143) 



24 



where 



Cf = ^^1^, Tf = Tnnj = ^Uf, = ^Nc - ^Tf (144) 



c 



and 



AP(J) = 2Tfi2x-l) 

AP(;) = 2iV,(^^^-2a; + l)+5(l-x)|, (145) 
while the LO transverse kernels are given by 

^-n':'=C.((^-2 4*(l-.)). (146) 
The unpolarized kernels should be compared with the Collins-Qiu form 

Pqq = Iqq-Si'^-X) dz'Jgq 

Jo 

= - dZ-fgg dZ-fgg^ 5(l " X) 



Pq9 '^19 

P99 ~ 7sij 



where 



Iqq = Cf(-^ l-X 

\1 ~ X 

Iqg = (2X - 1) 

\l — X X 



(147) 



(148) 



These kernels need a suitable regularization to be well defined. Below we will analize the 
implicit regularization underlying eq. (I147|) . One observation is however almost immedi- 
ate: the component Pgg is not of the form given by eq. p39|l . In general, therefore, in the 
singlet case, the generalization of eq. ()139j) is given by 



P{x)=Pi{x)-5{l-x) f p2{z)dz (149) 

Jo 

and a probabilistic interpretation is more complex compared to the non-singlet case and 
has been discussed in the original literature (Hj- 
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11 Convolutions and Master Form of the Singlet 



Distributions are folded with the kernels and the result rearranged in order to simplify the 
structure of the equations. Since in the previous literature this is done in a rather involute 
way 110] we provide here a simplificaton, from which the equivalence of the various forms 
of the kernel, in the various regularizations adopted, will be apparent. All we need is to 
apply ([7n|) and the evolution equations become 

dq „^ f dy yq{y) - xq{x) , „^ , s . . dy 3 



d\og{Q' 



f dy yq{y) - xq[x) , , . , . f dy ,^ . , ^ 3 

2Cf / h 2Cf log(l - X) q[x) - — [l + z) q[y) + -Cpqix) 

J y y - X y 2 



+-/ - [z' + (1 - zf) 9{y) 

J X y 



dg ^ r^dyl + {l-zy dyyf{y)-xf{x) 



ciiog(Q2) y z Jx y y-x 

^'dy fl _ ^^ , ^ , 



+2N, log(l - x)g{x) + 2iVc ^ ^ Q - 2 + ^(1 - z)^ g{y) + 



(150) 



with z = x/y. The same simplified form is obtained from the probabilistic version, having 
defined a suitable regularization of the edge point singularities in the integrals over the 
components 7/// in eq. ()148|) . The canonical expressions of the kernels ()145p . expressed in 
terms of "+" distributions, can also be rearranged to look like their equivalent probabilistic 
form by isolating the edge-point contributions hidden in their "+" distributions. We get 
the expressions 



(161) 



and 



the other expressions remaining invariant. 

A master form of the singlet (unpolarized) equation is obtained by a straightforward 
change of variable in the decreasing terms. We obtain 

dq f^^^ dy dy 

T = / —lqq{x/y)q{y)- / —Iqqiy/X)q{x) 

dr y Jo y 

J' = f —^99i^/y)-'^f[ iqgiy/^)9{x) 
X y 



1 dy 

2 l99{yl^)9{x)+ —lgq{x/y)q{y) (153) 
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with a suitable (unique) cutoff A needed to cast the equation in the form The 
(regulated) transition probabilities are then given by 



Wqg{x 

Wqq{y 

Wgg{X 

Wqqiy 
Wgq{y 



Wgq{X 



y) = iqq{x/y) 



X] 



X) 

x) 

y) 



lqq{y/ 



X] 



y) = lgg{x/y) 



9{y > x)d{y < 1 - A) 

y 

e{y < X - A)e{y > 0) 

X 

9{y > x)9{y < 1 - A) 

y 



nflqgiy/x) - \lggiy/x) \ ^ ^ ^^^^^ ^ 



Igqi^/y) 
0, 



e{y > x)e{y < 1 - A) 



(154) 



as one can easily deduct from the form of eq. p42jl . 



12 A Kramers- Moyal Expansion for the DGLAP Equa- 
tion 

A way to get rid of the integrals that appear in the master equation is the Kramers-Moyal 
(KM) expansion. These expansions (backward or forward) are sometimes useful in order 
to gain insight into the master equation itself, since they may provide a complementary 
view of the underlying dynamics. 

The expansion introduces differential operator of arbitrary order, due to the nonlocal 
structure of the equation. For the approximation to be useful, one has to stop the expan- 
sion after the first few orders and in many cases this turns out to be possibile. Examples 
of processes of this type are special Langevin processes and processes described by a 
Fokker-Planck operator. In these cases the probabilistic interpretation allows us to write 
down a fictitious lagrangean, a corresponding path integral and solve for the propagators 
using the Feynman-Kac fomula. For definitess we take the integral to cover all the real 
axis in the variable x' 



d_ 



q(x,T) 



dx' {w{x\x')q{x' , t) — w{x'\x)q{x, r)) dx' . 



(155) 



As we will see below, in the DGLAP case some modifications to the usual form of the 
KM expansion will appear. At this point we perform a KM expansion of the equation in 
the usual way. We make the substitutions in the master equation y ^ x — y m the first 
term and y x + y m the second term 



^g(x,r) 



dy {w{x\x — y)q{x — y,T) — w{x + y\x)q{x, r)) , (156) 



identically equal to 



dy {w{x + y — y'\x — y')q{x — y', r) — w{x + y'\x)q{x, r)) , (157) 



with y = y'. First and second term in the equation above differ by a shift (in —y') and 
can be related using a Taylor (or KM) expansion of the first term 
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-^q{x,T) = J dy22—j—^{w{x + y\x)q{x,T)) (158) 

where the n = term has canceled between the first and the second contribution coming 
from ()157|) . The result can be written in the form 

= E -^0^ («"(^)^(^' -)) (159) 

n=l 

where 

/OO 
dy{y - x)'^w{y\x). (160) 
-OO 

In the DGLAP case we need to amend the former derivation, due to the presence of 
boundaries (0 < x < 1) in the Bjorken variable x. For simplicity we will focus on the 
non-singlet case. We rewrite the master equation using the same change of variables used 
above 



^ q{x,T) = [ dyw{x\y)q{y,T)- I dyw{y\x)q{x,T) 



dr 



a{x) 

dyw{x + y\x) * q{x,T) + I dyw{x + y\x)q{x,T), (161) 



where the Moyal product gets simplified 



w{x + y\x) * q{x) = w{x + y\x)e ^q{x,T) (162) 



and a{x) = x — 1. The expansion is of the form 



j—x j-aix) ( — y)" 

/ dyw{x + y\x)q{x,T) — y^ / dy j — dx^ {w{x + y\x)q{x,T)) (163) 

Ja(x) ~i Jo ^■ 



which can be reduced to a differential equation of arbitrary order using simple manipula- 
tions. We recall that the Fokker-Planck approximation is obtained stopping the expansion 
at the second order 

d 1 

— g(x, r) = ao(x) - (ai(x)g(x)) + -dl (a2(x)g(x, r)) (164) 

with 

a„(x) = / dyy"'w{x + y,x) (165) 



being moments of the transition probability function w. Given the boundary conditions on 
the Bjorken variable x, even in the Fokker-Planck approximation, the Fokker-Planck ver- 
sion of the DGLAP equation is slightly more involved than Eq. ()164|) and the coefficients 
a„(x) need to be redefined. 
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13 The Fokker-Planck Approximation 



Having identified a probabilistic description of the DGLAP equation in terms of a master 
equation, it is then natural to try to investigate the role of the Fokker-Planck (FP) 
approximation to it. In the context of a random walk, an all-order derivative expansion 
of the master equation can be arrested to the first few terms either if the conditions 
of Pawula's theorem are satisfied -in which case the FP approximation turns out to be 
exact- or if the transition probabilities show an exponential decay above a certain distance 
allowed to the random walk. Since the DGLAP kernels show only an algebraic decay in x, 
and there isn't any explicit scale in the kernel themselves, the expansion is questionable. 
However, from a formal viewpoint, it is still allowed. With these caveats in mind we 
proceed to investigate the features of this expansion. 
We redefine 



r—x 

^0(2;) = / dyw{x + y\x)q{x,T) 

J a{x) 
faix) 

an{,x) = / dyy^'wix + y\x)q{x,T) 



/•a{x) 

dnix) = / dyy"-d^'' {w{x + y\x)q{x,T)) n = l,2, ... (166) 
^0 

For the first two terms (n = 1, 2) one can easily work out the relations 

di{x) = dxtti^x) — a{x)dxa{x)w{x + a{x)\x)q{x,T) 

d2ix) = dla2{x) — 2a{x)(dxa{x)Yw{x + a{x)\x)q{x, r) — a{xYdxa{x)dx {w{x + a{x)\x)q{x, r)) 
-a'^{x)dxa{x)dx {w{x + y\x)q{x, r)) \y=a{x) (167) 



Let's see what happens when we arrest the expansion (|163|) to the first 3 terms. The 
Fokker-Planck version of the equation is obtained by including in the approximation only 
dn with n = 0, 1, 2. 

The Fokker-Planck limit of the (non-singlet) equation is then given by 

d . . 1. 

— g(a;,r) = ao(x) + ai(a;) - -a2(x) (168) 

which we rewrite explicitly as 

d . . ^ /85 3 13 10 12 ^, , , 

A similar approach can be followed also for other cases, for which a probabilistic 
picture (a derivation of Collins-Qiu type) has not been established yet, such as for hi. We 
describe briefiy how to proceed in this case. 

First of all, we rewrite the evolution equation for the transversity in a suitable master 
form. This is possible since the subtraction terms can be written as integrals of a positive 
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function. A possibility is to choose the transition probabilities 

MAy] = —{—^-2]e{y>x)e{y<i) 

y X^-x/y J 



which reproduce the evolution equation for hi in master form 
dhi 



(170) 



dywi{x\y)hi{y,T) - dyw2{y\x)hi{x,T). (171) 



dr 

The Kramers-Moyal expansion is derived as before, with some slight modifications. The 
result is obtained introducing an intermediate cutoff which is removed at the end. In this 

case we get 

dr \ 6 Sx-^ x"^ X \ X J J 

^ (n 3 llx\ ^ , , , 

/3 1 5r^\ 
+Cf (2 - 3^ - 2^ + ^^^2) 

Notice that compared to the standard Fokker-Planck approximation, the boundary 
now generates a term on the left-hand-side of the equation proportional to q{x) which is 
absent in eq. (I164|) . This and higher order approximations to the DGLAP equation can 
be studied systematically both analytically and numerically and it is possible to assess 
the validity of the approximation 



13.1 Links to Fractional Diffusion 

A formal connection of the LO DGLAP equation to fractional diffusion can also be easily 
worked out |E2] in dimensional regularization. Since the literature on fractional calculus 
and anomalous diffusion is quite vast |2n|, we need just few essential definitions to make 
our brief discussion self-contained. The n-th primitive of a distribution function /(x) can 
be written as 

rf{x) = ^^^^ j\x - yr-^f{y)dy (173) 
which can be easily analytically continued for any real a > 

rf{x) = f f{y){x - yy-Hy (174) 

i ["J Jx 

thereby defining the /racizona/ integral of order a of the function f{x). We also recall that 
the fractional derivative of a function can be defined by a suitable analytic continuation 

where n — l<[3<n. In our case the role of (3 is taken by the parameter e of dimensional 
regularization. One can also formally define = J" to denote the corresponding 

fractional primitive. 
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Now, since the kernels contain "+" distributions and Dirac delta functions, is is con- 
venient to use the relation 



wj [i—w)+ \ L — w e 

valid in dimensional regularization, with e > 0, from which one obtains 

1 / 1 1 / 1 1 



6{l-w) + 0{e'), (176) 



2\l-wJ ^2{l-w) il-w)+' ^^^^^ 

The equation for the transverse spin distribution, for instance, can be easily recast in the 
form 

^ = -C^ {D^ + D-^) h{x) - 2CpD-' (M^) + Ic^hil) (178) 

with similar expressions for all the other non singlet equations. This equation is a possible 
starting point for the analysis of anomalous diffusion in the context of such equations. 
Equations of this type provide probability functions belonging to the class of stable dis- 
tributions, such as those describing Levy processes and continuous time random walks. 



14 Conclusions 

We have reviewed previous work of us on the determination of the solutions of renormal- 
ization group equations in QCD from x-space and on the KM expansion of the master 
formulation of these equations in analogy to the theory of stochastic processes. Parton 
distributions represent the first important example of nonlocal operators in a realistic 
field theory that have been studied extensively for almost three decades. Although we 
are unable to compute from first principles these fundamental field theory structures, the 
study of their identification and classification at leading twist and higher has allowed to 
break substantial new ground in field theory and QCD phenomenology. 
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